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A multifrequency method based on the Matched Multifilter 
for the detections of point sources in CMB maps 
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ABSTRACT 

In this work we deal with the problem of simultaneous multifrequency detection of 
extragalactic point sources in maps of the Cosmic Microwave Background. We apply a 
linear filtering technique that uses spatial information and the cross-power spectrum. 
To make this, we simulate realistic and non-realistic flat patches of the sky at two 
frequencies of Planck: 44 and 100 GHz. We filter to detect and estimate the point 
sources and compare this technique with the monofrequency matched filter in terms of 
completeness, reliability, flux and spectral index accuracy. The multifrequency method 
outperforms the matched filter at the two frequencies and in all the studied cases in 
the work. 

Key words: methods: data analysis - techniques: image processing - radio contin- 
uum: galaxies - cosmic microwave background - surveys 



1 INTRODUCTION 

Over the last few years, a big effort has been devoted 
to the problem of detecting point sources in Cosmic Mi- 
crowave Background (CMB) experiments. The main reason 
is that modern CMB experiments have reached resolution 
and sensitivity levels such that their capability to estimate 
the statistics of CMB fluctuations at high multipoles is no 
longer limited by instrumental noise but by Galactic and 
extragalactic foreground contamination. Among extragalac- 
tic contaminant s, point sources are the most relevant , both 
in temperature jToffolatti et al]|l99Sl: iDe Zotti et al.|[l999l : 
iHobson et all Il999l: |Pe Zotti et alJ l2005h and in polarisa - 



tion ijTucci et alJ 12004. 20051: lLopez-Caniego et~atl [20091. 
Moreover, they are one of the most difficult contaminants 
to deal with and, at least in the frequency range spanned by 
CMB experiments, one of the most poorly known. It is there- 
fore mandatory to detect the maximum possible number of 
extragalactic point sources (EPS) and to estimate their flux 
with the lowest possible error before any serious attempt to 
study the CMB anisotropies. 

But EPS are not only a contaminant to get rid of. The 
study of EPS at microwave frequencies is very interesting 
from the standpoint of extragalactic astronomy. The next 
generation of CMB experiments will allow one to obtain of 
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all-sky EPS catalogues that will fill in the existing observa- 
tional gap in our knowledge of the Universe in the frequency 
range from 20 to roughly 1000 GHz. We expect to derive 
source number counts and spectral indices, to study source 
variability and to discover rare objects such as inverted 
spectrum radio sources, extreme gigahertz peaked spectrum 
sources (GPS) and high-redshift dusty proto-spheroids (see 
for ex ample the Planck Bluebook, iThe Planck Collaboration! 
2006). The first fruits of these new era of CMB experiments 
are already here: th e Wilkinson Microwa ve Anisotropy Probe 
(WMAP) satellite ijBennett et al.ll2003l ) has made possible 
the obtaining of the first all-sky point source catalogues 
above ~ 0.8-1 Jy in the 23-94 G Hz ran ge of frequencies 

iii tem perature ^Bennett et all ll2003ll: iHinshaw et al.l 

tOOi); IChen fc Wrightl J2008T): TWright et al.l [12 0091 ; 



lLopez-Caniego et al.l (12 0071 



Massardi et al. 
120090. ~ 



(2009) and 



polarisation l|L6pez-Caniego et al.1 120091') . being complete 
the la st three ones. In a recent paper iGonzalez-Nuevo et al.l 
(2008) have shown how these catalogues can be used 
to study the statistics of point sources in tha t range of 
frequencies. The Planck mission (jTauberl lioolil l will allow 
us to extend these catalogues down to lower flux limits and 
up to 857 GHz. 

We have mentioned that the detection and estimation 
of the flux of EPS are a difficult task. The main reason for 
this is that the many different types of EPS that are dis- 
tributed in the sky form a very heterogeneous set of objects 
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that do not have a common spectral behaviour. While other 
foreground contaminants follow a specific emission law that 
is approximately well known (or can be inferred from ob- 
servations) and that varies relatively slow and continuously 
across the sky, each extragalactic point source has an emis- 
sion law, that, in principle, can be totally different to any 
other and independent from them. Prom the point of view of 
statistical signal processing, the problem of detecting EPS 
is a case of sorely under-determined component separation 
problem where the number M of components is much larger 
than the number iV of frequency channels. 

Let us consider the case of another contaminant that 
is of relevance in CMB images and to which the EPS 
proble m has some similarities: the S uny a e v- Zel ' do vich ( S Z ) 
effect ijSunvaev fc Zeldovichl Il97d . Il972l : [ Carls trom et al.l 
2002). As in the case of EPS, the SZ contamination oc- 
curs in small and compact regions of the sky and it af- 
fects the high-^ region of the fluctuations angular power 
spectrum. But in the case of the thermal SZ effect the 
frequency dependence is very well known (ignoring rela- 
tivists corrections). This has made possible the develop- 
ment of a number of powerful detection techniques that 
are specifically suited t o the problem of SZ effect de- 
tection in CMB images llDiego et al.l 120021: iHerranz et al.l 
002al: iHobson fc McLachlanl 12003 riPierpaoli et all | 200a 
Herranz et all 120051: IVale fc Whitell200d: bires et al] bood : 
Melin et al.ll200d : iBartlett et al.ll20Qgl : iLeach et al.ll2008h or 



to use more generic diffuse component se paration techniques 
(see for example, among many others "iMaino et al. 20021: 
Stolvarov et al.ll2002l : lDelabrouille et~aTll2003| : lEriksen et al] 



Diva 

2004 : bedim et al.l [20051 : iBonaldi et alj|200d ) to obtain SZ 



maps and catalogues. 

Unfortunately, these approaches are not directly appli- 
cable to EPS detection. The most commonly used approach 
to this problem consists on working separately in each chan- 
nel. The key idea is to take advantage of the fact that all the 
EPS have the same shape (basically, that of t he beam). In 
the fi el d of C M B images, wavelet tech n iques llVielva et al 



2001 



2006 
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120031: iGonzalez Nuevo et all l200d : ISanz et al 
Lopez-Canieg o et al] 120071). matched filters (MF , 

iBarreiro et"aTI 120031 : 
fil- 



Lopez-Caniego et al]_ 200d) and oth e r related line ar 



tering techniques |Sanz et al] 2001 ;_ Chiang et al] 



Herranz et al] 002cj ; lLopez-Caniego et al] 20041 . 



2002 



2005 



Lopez-Caniego et al] 120051 ) have proved to be useful. All 



these techniques rely on the prior knowledge that the 
sources have a distinctive spatial behaviour and this fact 
is used to design some bandpass filter to enhance them 
with respect to the noise. Detection can be further im- 
proved by including prior information about the sources, 
i.e. some knowledge about th eir flux distribution, in the 
frame of a Bayesian formalism ijHobson fc McLachlanll2003l : 
ICarvalho et~afll2009l ). 

With the previous single frequency methods and for the 
case of the Planck mission we expect to reach detection limit 
fluxes that range from a few hundred mjy to several Jy, de- 
pending on the frequency channel, and to obtain catalogues 
that wil l contain from several hundreds to a few tho usand 
objects ijLopez-Caniego et alJbood : ILeach et alj|2008l ). This 
is satisfactory, but we feel that we could do better if we were 
able to use multi-wavelength information in some way. 

For now, multi-wavelength detection of EPS in 



CMB images remains a largely unexplored field. In re- 
cent years , some attem p ts ha v e been done i n this 
direction liNaselskv et~a71 120021 ; IChen fc Wrightl 120081 : 
IWright et all 120091 ). More recently. IHerranz fc Sanzl (2008) 
have introduced the technique of 'matched matrix filters' 
(MTXF) as the first fully multi-frequency, non-parametric, 
linear filtering technique that is able to find EPS and to 
do unbiased estimations of their fluxes thanks to the dis- 
tinctive spatial behaviour of the sources, while at the same 
time does incorporate some multi-wavelength information, 
without assuming any specifi c spectral behaviour for the 
sources. IHerranz et al] l|2009l ) have applied the MTXF to 
realistic simulations of the Planck radio channels, showing 
that it is possible to practically double the number of detec- 
tions, for a fixed reliability level, for some of the channels 
with respect to the single frequency matched filter approach. 

MTXF use multi-wavelength information in such a way 
that it is not necessary to make any assumption about the 
spectral behaviour of the sources. In fact, in that formalism 
their spectral behaviour is entirely irrelevant. All the multi- 
wavelength considerations concern only to the nois^B and its 
correlations. In this sense, MTXF deal only with half of the 
problem. This has its advantages in terms of robustness and 
reliability, but one could wish to have a technique that uses 
multi-wavelength information in the modelling of both the 
signal (EPS) and the noise. But, as we mentioned before, 
the spectral behaviour of the EPS is not known. 

In this paper we will show that even if the spectral be- 
haviour of the EPS is unknown a priori, it is still possible to 
determine it directly from the data by means of an adaptive 
filtering scheme that incorporates multi-frequency informa- 
tion not only through the noise correlations among channels, 
but also about the sources themselves. 

The problem is, in more than one sense, similar to the 
problem of detecting SZ clusters (that is the reason we dis- 
cussed that case a few paragraphs above). In the SZ case 
the spectral behaviour of the sources is known, but not 
their size. A way to deal with this is to introduce the scale 
of the source as a free parameter (for example, the clus- 
ter core radius r c ) in the design of a 'matched multifilter' 
(MMF) and to optimise the value of this parameter for 
each source so that a maximum sign al to noise ratio is o b- 
tained after filtering (see the details in[Herranz et al.ll002al lb1. 
120051 : ISchafer et al"Tl200d : iMelin et al.ll200d ). As the prob- 
lem depends on the optimisation of one single parameter, 
the method is easy to implement in codes that are relatively 
fast. 

In this paper we introduce a modification of the MMF 
technique in which the mixing coefficients of the frequency 
dependence vector of the sources are considered as free pa- 
rameters to be optimised. As we will show, if the number 
of frequency channels is N and we choose wisely a fiducial 
frequency of reference, the number of free parameters to op- 
timise is N— 1. For simplicity, throughout this paper we will 
use as an example the case of two channels, N = 2, but the 
method is valid for any number of channels N > 1. 

The structure of this paper is as follows. In section[2]we 
will summarise the formulae of the MMF and we will intro- 



1 Here the term 'noise' refers to all the components except for 
the EPS themselves, including CMB and Galactic foregrounds. 
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duce the modification that allows us to use it for detecting 
EPS in arbitrarily chosen pairs of frequency channels. We 
will show that the filter can be in this case factorised in 
such a way that a particularly fast implementation of the 
method can be achieved. In section [3] we will describe the 
simulations that we use to test the method. The results of 
the exercise will be described in section [4] Finally, in sec- 
tion [5] we will draw some conclusions. 



2 METHOD 

2.1 The single frequency approach 

Let us assume a set of images corresponding to the same 
area of the sky observed simultaneously at N different fre- 
quencies: 



y„(x) = f„s v (x.) +n„(x), 



(1) 



where v = 1,...,N. At each frequency v, y v is the total 
signal in the pixel x and s v represents the contribution of 
the point source to the total signal y,y, for simplicity let 
us assume there is only one point source centered at the 
origin of the image; /„ is the frequency dependence of the 
point source; and n u is the background or generalized noise 
(containing not only the instrumental noise, but also the 
contributions of the rest of components). 

The intrinsic angular size of the point sources is smaller 
than the angular resolution of the detector. At each observ- 
ing frequency, each source is convolved with the correspond- 
ing antenna beam. For simplicity we will assume the antenna 
beam can be well described by a symmetric 2D Gaussian 
function. Then we can write 



s v (x) = At v {x), 



(2) 



where x = |x| (since we are considering symmetric beams), 
A is the amplitude of the source and r is the spatial template 
or profile. The background n„(x) is modelled as a homoge- 
neous and isotropic random field with average value equal 
to zero and power spectrum P v defined by 



(n„(q)n*(q')) = P„<5o(q - 



(3) 



where n„(q) is the Fourier transform of n„(x) and 8 2 D is the 
2D Dirac distribution. 

In the single frequency approach each channel is pro- 
cessed separately and independently from the other fre- 
quency channels. This approach is robust in the sense that 
it is not necessary to assume anything about the spectral 
behaviour of the sources. The main drawback of the single 
frequency approach, however, is that one misses the poten- 
tial noise reduction that could be obtained with a wise use 
of the information present at the other frequencies. 

The standard single frequency point source detec- 
tion method in the literature is based on the matched 
filter llTegmark fc de Oliyeira-Costal Il998l : iBarreiro et al.l 
120031 : lLopez-Caniego et al.ll2006l ). The matched filter is the 
optimal linear detector for a single map in the sense that 
it gives the maximum signal to noise amplification. The 
matched filter can be expressed in Fourier space in the fol- 
lowing way: 



</>mf(<?) 



r(q) 

a p( q y 



dq 



P(q)' 



(4) 



Here a is a normalisation factor that preserves the source 
amplitude after filtering. 



2.2 The Matched Multifilter 

In the multi-frequency approach we take into account the 
statistical correlation of the noise between different fre- 
quency channels and the frequency dependence of the 
sources. Now let us model the background n„(x) as a ho- 
mogeneous and isotropic random field with average value 
equal to zero and crosspower spectrum P vi v 2 defined by: 



K(q)n*,(q')) =-FV 2 <5!>(q-q') 



(5) 



where n„ (q) is the Fourier transform of n u (x) and 8% is the 
2D Dirac distribution. Let us define a set of N linear filters 
ip v that are applied to the data 



io„(b) 



dx j/„(x)y>„(x; b) 



dq e * q l V(q)?/v(<7). (6) 



Here b defines a translation. The right part of equation (H} 
shows the filtering in Fourier space, where yv{(\) and i)iv(q) 
are the Fourier transforms of y v {x) and ^„(x), respectively. 
The quantity w v {bi) is the the filtered map v at the position 
b. The total filtered map is the sum 



<;(b)=5>„(b)- 



(7) 



Therefore, the total filtered field is the result of two steps: 
a) filtering and b) fusion. During the first step each map y v 
is filtered with a linear filter ip v ; during the second step the 
resulting filtered maps w v are combined so that the signal s 
is boosted while the noise tends to cancel out. Note that the 
fusion in eq. (J7)) is completely general, since any summation 
coefficients different than one can be absorbed in the defini- 
tion of the filters ij) v . Then the problem consists in how to 
find the filters %j) v so that the total filtered field is optimal 
for the detection of point sources. 

The total filtered field w is optimal for the detection of 
the sources if 

(i) iti(O) is an unbiased estimator of the amplitude of the 
source, so (w(0)) = A; 

(ii) the variance of u)(b) is minimum, that is, it is an 
efficient estimator of the amplitude of the source. 

If the profiles t v and the frequency dependence /„ are known 
and if the crosspower spectrum is known or can be estimated 
from the data, the solution to t he problem is already known 
the matched multifilter (MMF, iHerranz et 



3 already . 
aT1 l002ah : 



/ 



dq F'P _1 F, 



(8) 



where ^f(q) is the column vector ^f(q) = [i/v(g)], F is the 
column vector F = [fuT v ] and P 1 is the inverse matrix 
of the cross-power spectrum P. Finally, we can obtain the 
variance of the total filtered field, given by the following 
expression: 



/ 



dq* P* = a 



(9) 
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2.3 MMF with unknown source frequency 
dependence 

As it was previously discussed, the problem is that the fre- 
quency dependence f v of the sources is not known a pri- 
ori. Then, the possibilities are either a) to admit defeat, 
returning to the single frequency approach, b) to devise a 
filtering method that does not use the frequency depen- 
dence of the sources altogether or c) to model somehow 
the unknown frequency dependence in the framework of 
some opt imisation scheme . The second approach was ex- 
plored in iHerranz fc Sand <|2008h : IHerranz et al.l <(200Sh . In 
this work we will study the third approach to the problem. 

Before addressing this problem, it will be useful to 
rewrite eq.© in a slightly different way. Let us write the 
vector F = [/i/7v] in matrix form as 

F = T{q)f{v), (10) 

with T a diagonal matrix T(q) = diag[ri (q), . . . , Tff(g)] and 
f = [/„] the vector of frequency dependences. Note that all 
the dependence in q is included in the matrix T; this fact 
will be useful later. 

Now imagine that f describes the true (unknown) fre- 
quency dependence of the sources and that g = [g v \, v = 
l,...,N is a new vector, of equal size as f but whose ele- 
ments can take any possible value. We can define the MMF 
for vector g 

* g (g) = a g P^Tg, 

"g 1 = J dq g'TP-'Tg = g'Hg, (11) 

where H = / dqTP _1 T and we have used the facts that 
T = T and that vector g does not depend on q and can 
therefore go out of the integral. When applied to a set of 
images where there is present a point source with true am- 
plitude A and true frequency dependence f , the filters 
will lead to an estimation of the amplitude 

A s = w s (0) = a s A g'Hf. (12) 

Note that if g 7^ f , then A s 7^ A. On the other hand, the 
variance of the filtered field would be, in analogy with eq. 
([9]), o-g = a g . Let us finally define the signal to noise ratio 
of the source in the total filtered map as 

SNR S = ^2.. (13) 

Then we can ask what is the vector g that maximizes the 
signal to noise ratio SNR S . Intuition alone indicates that 
SNR S is maximum if and only if g = f . This can be formally 
proved with little effort by taking variations of g. 

Then the problem can be solved via a maximisation al- 
gorithm. In the case of a non blind search, where the position 
of a given point source is known, one can focus on that point 
source and iteratively try values of the elements of g until a 
maximum signal to noise is reached. In the case of a blind 
search, the situation is a little bit more difficult because in 
a given image there may be many different objects S; with 
a different solution gi. A way to proceed is to filter many 
times the image, using each time a different set of values 
of the elements of g so that the appropriate range of fre- 
quency dependences is sufficiently well sampled, and then 
to proceed counting one by one all the possible detections 



and associating to each one the values of g that maximize 
the signal to noise ratio of that source in particular. 

We would like to remark that this situation is very sim- 
ilar to the case of the detection of galaxy clusters wi th un- 
known angular size described in IHerranz et al ] (|002al fbh. In 
that case the frequency dependence f was known but the 
size of the clusters (their source profile) was not. The clus- 
ter profile can be parametrised as a modified beta profile 
with a fr ee scale parameter r e (typically, the cluster core ra- 
dius). In IHerranz et al.l l|002al lbh it was shown that the true 
scale of the clusters can be determined by maximizing the 
signal to noise ratio of the detected clusters as a function 
scale r c of the filter. 

In our case, factorisation ifTOll leads to equations (fTTT) 
and ll"12ll : this is very convenient for implementation of the 
MMF when many filtering steps are necessary. The most 
time-consuming part of the filter is the calculation of matri- 
ces P and T because they must be calculated for all values 
of q. In the case of clusters with unknown size T had to 
be calculated for every value of r c . However, in the case we 
are considering in this paper the only quantity that varies 
during the maximisation process are the elements of vector 
g. This allows us to compute the integrals of matrix H only 
once for each set of images. As a result, applying the MMF 
to large numbers of point sources with unknown frequency 
dependence is, in general, much faster than applying the 
MMF to the same number of clusters with unknown source 
profile. 

The main difference is that while in the case of galaxy 
clusters it was necessary to maximize with respect to only 
one single parameter (the core radius), in the case of the 
unknown frequency dependence it is necessary to maximize 
with respect to the N components of vector g. This proce- 
dure can require a very large number of computations if N 
is big. Although we have just seen that each free parame- 
ter of the frequency dependence can be mapped much faster 
than each free parameter of the source profile, we are still 
interested in reducing the number of computations as much 
as possible. 



2.4 Number of degrees of freedom of vector g 

For N images, vector g = [51,... , gjv] has N degrees of free- 
dom. This makes the optimisation procedure more complex 
and computationally expensive. The situation can be light- 
ened if we choose one of the frequencies under consideration 
as our fiducial frequency of reference. Let us choose for ex- 
ample a concrete frequency j £ {1, . . . , N} to be our fiducial 
frequency of reference, then 

( yj (0)) = Af jTj (0) = A (14) 

and therefore, since the profile Tj is normalised to unity, fj 
must be equal to one. Therefore, we must look for vectors 
g = gi, . . . , gj-i, 1, gj+i, ■ ■ ■ , N and the number of indepen- 
dent degrees of freedom is N — 1. If the number of channels 
is N = 2, there is only one degree of freedom for the op- 
timisation problem. In the next sections of this paper we 
will consider, for simplicity, the case of two channels, but we 
would like to remark that the extension to N > 2 frequencies 
is straightforward. 
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2.5 Optional parametrisation of vector g 

Another way to reduce the number of degrees of freedom 
of vector g is to find a suitable parametrisation for it. For 
example, the power law relationship 

7H = /o(£)"\ (15) 

where I{v) is the flux at frequency v, vq is a frequency of 
reference, Jo is the flux at that frequency of reference and 
7 is the spectral index, is widely used in the literature. If 
eq. (1 1 5 j) . then the reference flux Jo can easily be related to 
the reference amplitude A of the sources and the number of 
degrees of freedom is just one, the spectral index 7. 

However, the parametrisation of vector g has its own 
risks. For example, it is known that eq. lllSIl is valid only 
as an approximation for any given frequency interval and 
that its validity decreases as the size of the interval grows. 
If we choose to follow the parametric approach we know for 
sure that the results will be less and less accurate as the 
number N of channels grows, especially if the separation 
between frequency channels is large. On the other hand, eq. 
lfl~5]l is alway exact if = 2. Therefore we can safely use 
the parametrisation lfl~5]l for N = 2 channels, without loss 
of generality. Since the number of degrees of freedom is one 
either if lfl~5]l is used or not, the use of the parametrisation 
is irrelevant in this case. However, we may be interested in 
using it for historical, didactic and practical motivations. For 
example, eq. lfl~5]l is useful to express the physical properties 
of the sources in terms of their (steep, flat, inverted, etc.) 
spectral index. 

We would like to remark that frequency dependence 
parametrisation, in the form of eq. lfl~5]l or any else other way, 
may or not be useful in some cases, but it is not essential 
at all for the method we propose in this paper. 



3 SIMULATIONS 

In order to illustrate the MMF method described above and 
to compare the MMF multi-frequency approach with the sin- 
gle frequency approach, we have performed a set of basic, yet 
realistic, simulations. We consider the case of two frequency 
channels (N = 2). Generalisation to more frequency chan- 
nels is possible, as discussed in section |2~3"1 but we choose to 
keep things simple in this paper. 

For this exam ple we take the case of the Planck mis- 
sion l|Taubenl2005l ). We will consider the 44 GHz and 100 
GHz Planck channels. The choice of the pair channels is not 
essential: any other pair of channels would have served the 
same for this exercise. This particular choice allows to study 
the case of radio sources in two not adjoining channels with 
different instrumental settings: the 44 GHz channel belongs 
to the Low Frequency Instrument of Planck and the 100 
GHz channel belongs to the High Frequency Instrument. 

For the si mulations we have use d the Planck Sky 
Mode0 (PSM, iDelabrouille et al.l liooil . in preparation), a 
flexible software package developed by Planck WG2 for 



2 http://www.apc.univ-paris7.fr/APC_CS/Recherche/Adamis/ 
PSM/psky-en.php 



making predictions, simulations and constrained realisa- 
tions of the microw ave sky. The simul ated data used here 
are the same as in iLeach et al.l l|2008l ). where the charac- 
teristics of the simulations are explained in more detail. 
Maps are expressed in (AT/T), thermodynamic units. Sim- 
ulations include all the relevant astrophysical components: 
the CMB sky is based on a Gaussian realisation assum- 
ing the WMAP best-fit Ct at higher multipoles; Galactic 
emission is described by a three component model of the 
interstellar medium comprising free-free, synchrotron and 
du st emissions. Fr e e-free emission is based on the model 
of iDickinson et al.l l|2003l ) assuming an electronic temper- 
ature of 7000 K. The spatial structure of the emission is es- 
timated using a Ha template corrected for dust extinction. 
Synchrotron emission is bas e d on a n extrapolation of the 408 
MHz map of lHaslam et~aT1 jl982h from which an estimate 
of the free-free emission was removed. A limitation of this 
approach is that this synchrotron model also contains any 
dust anomalous emission seen by WMAP at 23 GHz. The 
thermal e mission from interstellar dust is estimated using 
model 7 of iFinkbeiner et alJ l|l999n . 

For the purely descriptive purposes of this example, 
we take eight different regions of the sky located at in- 
termediate Galactic latitude (four of them uniformly dis- 
tributed across the 40° North Galactic latitude parallel and 
four of them distributed in the same way 40° South of the 
Galactic plane). For each region we select a 512 x 512 pixel 
square patch (at 44 and 100 GHz). Pixel size is 1.72 ar- 
cmin for the two frequencies. Therefore, each patch covers 
an area of 14.656 square degrees of the sky. When both 
patches have been selected, we add simulated extragalac- 
tic point sources with a spectral behaviour described by eq. 
I|15p . We take as frequency of reference vq = 44 GHz. Note 
that eq. Illoll is expressed in flux units and the maps are in 
(AT/T) t h: we make the appropriate unit conversion before 
adding the sources. The antenna beam is also taken into 
account: the full width at half maximum is FWHM=24 ar- 
cmin for the 44 GHz channel and FWHM=9.5 arcmin at 
100 GHz. Finally, after doing that, we have added to each 
patch unifo rm white noise with the nominal levels specified 
for Planck ijThe Planck Collaboration! I2OO6) and this pixel 
size. 

We are interested in comparing the performance of the 
multi-frequency approach with that of the single-frequency 
matched filter. In particular, we expect to be able to de- 
tect fainter so urces with the MMF than with the MF. From 
recent works ijLopez-Caniego et al.ll200^ : ILeach e t al. 2008) 
we know that in this kind of Planck simulations, the MF 
can detect sources down to fluxes ~ 0.3 Jy (the particular 
value depends on the channel and the region of the sky). 
Here we will simulate sources in the interval [0.1, 1.0] Jy 
plus a few cases, that will be described below, where even 
lower fluxes are necessa ry. Regarding the spec t ral in dex of 
the sources, according to lGonzalez-Nuevo et al.ll|2008h . most 
radio galaxies observed by WMAP at fluxes ~ 1 Jy show 
spectral indices that lie in the range ( — 1.0, 1.4). 

We sample the interesting intervals of flux and spec- 
tral index by simulating sources with fluxes at 44 GHz Jo = 
{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0} Jy and spectral in- 
dices 7 = {-1.0,-0.7,-0.4,-0.1,0.2,0.5,0.8,1.1,1.4}. For 
each pair of values (Jo, 7) we have simulated 100 point 
sources. The point sources are randomly distributed in the 
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maps (obviously, the same source is placed in the same pixel 
in both frequencies), with only one constraint: it is forbid- 
den to place a source closer than FWHM44/2 pixels from 
any other. In this way we avoid source overlapping. Image 
borders are also avoided. For each set of 100 sources we pro- 
ceed in the following way: we randomly choose one among 
the eight patches we have and place 10 sources in it. Then we 
randomly choose other patch (allowing repetition) and place 
the next 10 sources, and so on. In total, we have simulated 
9000 sources for this exercise (the additional simulations at 
fluxes below 0.1 Jy are not included). 



4 RESULTS AND DISCUSSION 

In order to compare the matched filter and the matched mul- 
tifilter, we use the same maps with both methods. It means 
that not only the maps, but the sources are identical for the 
two filters (their intrinsic fluxes and positions). Each simu- 
lation is filtered separately with the matched filter (UJ) and 
the matched multifilter |[8|). In order to have a better esti- 
mation of the power spectra, avoiding as much as possible 
aliasing effects, we implement a power spectrum estimator 
that uses the 2D Hann window ij Jiang et al]|2002T ). 

We will compare performance of the two methods in 
terms of the following aspects: spectral index estimation, 
source detection and flux estimation. 

4.1 Source detection 

Direct comparison of source detection between the MF and 
the MMF is not an obvious task because for N — 2 input 
images the MF produces two filtered images, whereas the 
MMF produces only one combined filtered map. Whereas 
the meaning of 'detection' in the latter case is straightfor- 
ward (once a detection criterion is chosen, a source is de- 
tected or not), in the former the situation is not so clear. 
Imagine we decide to apply the same detection criterion to 
the two MF filtered images (which is not an obvious option), 
then for a given source we can have three different outcomes 
of the detection: 

• The source may be detected in both maps. 

• The source may be detected in only one of the maps. 

• The source may be detected in none of the maps. 

In the first two cases we can obtain point source catalogues 
(with different number of objects, in principle, for the two 
frequencies), but only in the first case are we able to estimate 
the flux at the two frequencies and therefore the spectral 
index. In the case of the MMF, if the source is detected we 
automatically know the spectral index and we can use eq. 
111. j|| to give the fluxes at the two frequencies. 

Therefore, in the following we will distinguish two dif- 
ferent cases when we speak about detections with the MF. 
On the one hand, the intersection of the detections in the 
two channels gives us the objects that can be used for study- 
ing the spectral index distribution; on the other hand, the 
union of the two sets gives us the total number of objects 
that can be detected in, at least, one of the channels. For 
the MMF, both sets are the same by definition. 

Regarding the detection criterion, for simplicity we 
will apply the same criterion to all the filtered maps: the 



widespread 5a threshold. Note that the 5cr threshold corre- 
sponds to different flux values for different filters. However, 
in this paper we will follow the standard 5a criterion for 
simplicity. 

Figure [T] shows the real sources (in %) that we detect 
above a 5a level detection whose intrinsic fluxes (values in- 
troduced by us in the simulations) in the reference frequency 
(Jo) are the corresponding values in the horizontal axis. Ta- 
bled] shows the flux at which we are able to detect, at least, 
the 95% of the sources. We can observe several interesting 
aspects. The first one is the fact that the matched multifilter 
improves the level of detection with respect to the matched 
filter level for all the values of 7 we have inserted. 

The second one is a natural selection effect: we detect 
more flat/inverted sources (7 ~ 0/ negative values of 7) at 
low fluxes than steep ones (positive values of 7). Figure [T] 
and Table [T] show us that the level of detections is higher 
for negative values of 7 than for positive values. Keeping in 
mind that the reference frequency vq is equal to 44 GHz, 
and according to eq. (j 1 5 1> . it can be seen that for 7 > 0, the 
simulated sources satisfy this condition: J100 < J«- In this 
case, the sources appear less bright at 100 GHz. In these 
conditions it is quite difficult to detect sources at 100 GHz. 
Therefore, it means that we are not able to give the spectral 
indices of these point sources by means of the matched filter 
method when 7 is strongly positive (for instance, 7 > 1). 
Obviously, the smaller 7 is, the better the detection is with 
the matched filter at 100 GHz. Therefore, we add from 7 = 
—0.1 two additional bins at Jo = 25, 50 mjy. 

Another aspect we have to remark is the similar aspect 
of the detection curve for the matched filter at 44 GHz for 
all the 7 values (see Figure [l}. The reason of this similar- 
ity is that the reference frequency is 44 GHz, and the maps 
we have simulated at this frequency are the same, indepen- 
dently of 7, with only one exception: the position of the 
sources. This means that statistically are equivalent (with 
the inherent fluctuations due to the variation in the posi- 
tions of the sources). We find a similar number of detections 
by means of the matched filter at 44 GHz, independently of 
7. For this reason it is difficult to characterize the spectral 
behaviour of the sources for Jo < 0.5 Jy, because we do not 
have a high percentage of detected sources at 44 GHz below 
that value of 7o. 

Additionally, we observe that the matched multifilter is 
capable to detect sources whose Jo < 0.1 Jy for 7 < —0.1. It 
is interesting to compare this with the matched filter, that 
does not detect sources below 0.1 Jy in the conditions of 
this work. This is a result that we obtain with the method 
presented here. It allows us to detect point sources whose Io 
is too low to be detected with the traditional matched filter. 

To summarise, we can say that the MMF improves the 
detection level. Specially remarkable are the cases where the 
sources are near to be flat (central row of Figure H). At 100 
GHz, the MF recovers the 100% of the sources for Jo ~ 
0.4 — 0.6 Jy. Meanwhile, the MMF reaches this level for 
Jo ~ 0.1 Jy. This particular case is really interesting in the 
sense that most of the sources have this spectral behaviour. 

4.2 Spectral index estimation 

As it was mentioned before, one of the quantities we want 
to obtain is the spectral index of the sources (eq. I15|l . As we 
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Figure 1. Number of detections against the input value of Io for different values of the spectral index 7. MF44 represents the sources 
detected with the matched filter at 44 GHz. MFIOO the same but at 100 GHz. MFi is the intersection of MF44 and MFIOO, and MFu is 
the union of MF44 and MFIOO. 



can see in that equation, when we know this spectral index 
and the flux at the reference frequency, we are able to give 
an estimation of the flux at other channels. In Figure [2] we 
see how we recover the spectral index by means of the MMF, 
and compare these results with the obtained by the matched 
filter. In general, we can observe that we can reach lower 
values of Io with the MMF. Also, at higher values of Io than 



> 0.4 Jy, we are able to give, with a good degree of precision, 
an estimation of 7 by means of the multifrequency method. 
In general, the error bars at these values of Io are quite 
smaller than the bars of the matched filter, so we recover 
with more accuracy the spectral index and less uncertainty. 

Other aspect is that the error bars increases when Io 
is smaller. It seems logical, because we have fainter sources 
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Table 1. Fluxes in the reference frequency {Iq) for which we detect, at least, the 95% of the sources for the different filtering methods. 
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Figure 2. Values of 7 recovered by means of the MMF (asterisks) and the MF (circles). The line indicates the ideal recovering. The 
circles corresponding to the MF are slightly displaced in the horizontal axis in order to distinguish the results. 
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and a smaller number of detections (see Figure [TJ. Then, 
at Jo = 0.1 Jy, we can see that the estimation of 7 is not 
as good as we wish, because it has a great uncertainty. The 
main reason is that the signal to noise ratio is close to the 
threshold level we have imposed. 

With the matched filter we estimate correctly the spec- 
tral index for To > 0.7 Jy at 7 ^ 0.8. At higher values of 
7 we find the same problem that we have mentioned be- 
fore: there are not so many detections at 100 GHz below 
0.7 Jy (Figure HJ. Since the detected sources are close to 
the noise level, the fluxes recovered present an overestima- 
tion with respect to t he input value due to the Eddington 
bias l|Eddingtonlll913h . an effect produced close to the noise 
level by the overestimation due to the fluctuations because 
of the noise in the positions where the source is located. 
As we said before, and seeing the Figure [TJ we detect more 
sources at 100 GHz than at 44 GHz for values of the spectral 
index smaller than 0.8. 

Finally, we can observe an interesting aspect of the 
matched filter. When we do not have the sufficient detections 
in, at least, one channel (the sources detected are below the 
~ 40% of the total number of sources), the estimation of 
the spectral index is not good. In all cases we see an over- 
estimation of 7, except for 7 = 1.4, 1.1, where we have an 
underestimation. This is due to the fact that at these values 
of 7 the sources at 100 GHz are fainter and the number of 
detections at this channel is really small. Then, because of 
the Eddington bias, the flux at this frequency is overesti- 
mated, and consequently, the value of 7 is underestimated. 
The Eddington bias explains as well the overestimation of 7 
in the other cases. The only difference is that now, it is at 
44 GHz where we have a smaller number of detections. If we 
also see the Figure [TJ we observe that for values of To < 0.6 
Jy, we are pretty close to the noise level. It means that the 
noise fluctuations in the maps produce an overestimation in 
the flux at 44 GHz (7o) and, in this case, an overestima- 
tion in 7 too. Summarising, for 7 = 1.4, 1.1, the Eddington 
bias appears at 100 GHz (underestimation of the spectral 
index). For the rest of values of 7, this bias appears at 44 
GHz (overestimation of the spectral index). 

4.3 Flux estimation 

The other quantity that we want to recover is the flux of the 
source at the reference frequency Vo- As we said in section 
14.2 j and according to eq. [15] when we obtain the flux at uq, 
we are able to estimate the flux at 100 GHz. 

Figure [3] shows the recovered flux at the reference fre- 
quency (44 GHz) for a given value of the spectral index. The 
error bars recovered with the matched filter are, in general, 
larger than the ones we obtain with the matched multifil- 
ter. It is particularly notorious at small values of To, where 
the recovered values of the flux have a good agreement with 
respect to the input values, with small error bars. 

We observe in Figure [3] the evolution of the recovered 
flux. In general, for all the values of 7 that we have studied, 
the matched multifilter is a suitable and effective tool to 
estimate the To of the sources. For the matched filter, we 
observe a good determination of To for input values above 
0.7 Jy. For smaller values, To has a higher value than its real 
one. That is due to the Eddington bias at 44 GHz. At this 
frequency, in Figure[T]we observe that for values smaller than 



0.6 Jy, we only detect a ~ 40% of the total sources. That 
means that many of these sources are close to this noise level. 
And for the correct estimation of To and the spectral index, 
we need a good detection of the sources at the two channels. 
For low values of To the number of detected objects is small 
and we have few statistics. 

In this section we have discussed about the flux at the 
reference frequency (44 GHz in our case). But we have used 
the matched multifilter with two different frequencies. For 
this reason it is necessary to say something related to the 
second channel at 100 GHz. It is important to obtain the 
values that we recover at 100 GHz with the matched mul- 
tifilter because the corresponding errors bars of To and 7 
could propagate additional errors in the flux estimation at 
100 GHz (see eq. lfl~5jl). After extrapolating the results, we 
obtain the flux at 100 GHz, and compare it with the results 
obtained with the matched filter. As we saw at 44 GHz, the 
improvement with the matched multifilter is clear respect to 
the matched filter. 

4.4 Reliability 

For academic purposes, in the previous sections, we have 
produced the simulations introducing 100 sources for each 
of the pair values of intensity and spectral index (see sec- 
tion [3}, that simplifies the comparison between both fil- 
ters for all the cases under study. On the contrary, it is 
well known that the number of sources pe r flux interval, 
the source number counts, is not constant ijDe Zotti et al.l 
120051 : iGonzalez-Nuevo et al.l lipoH) nor the spectral index 
distribution llSadler et al]|2008l : lGonzalez-Nuevo et al.ll2008l : 
iMassardi et ail 20091 In order to study the performance of 
the new method under more realistic simulations we pro- 
duced a new set of simulations (100) with the following char- 
acteristics: 

• We used as a background the same eight regions de- 
scribed in the previous sections. 

• The sources were simulated with an almo st Poissonian 
distribution (see IGonzalez-Nuevo et all l|2005h for more de- 
tails about the method) at 44GHz, with fluxes that foll ow 
the source number counts model of iDe Zotti et all ([2005) . 

• The fluxes at 100GH z were estimated assum i ng ran dom 
spectral indices from the IGonzalez-Nuevo et al.l l|2008l ) dis- 
tribution. 

• The point source maps were filtered with the same reso- 
lution as the background maps and randomly added to them. 

There is also another interesting quantity commonly used in 
the study of the performance of a source detector: the num- 
ber of spurious sources. Spurious sources are fluctuations of 
the background that satisfied the criteria of the detection 
method and therefore are considered as a detected sources. 
It is clear that the best method will be the one that has the 
best detections vs. spurious ratio. Therefore, this time we 
use more realistic simulations and count the spurious and 
real sources: the maps are filtered using the MMF and the 
MF at both frequencies, we estimate the position and inten- 
sity of the sources above 3<r level and, by comparing with 
the input source simulations, we count the number of real 
and spurious sources that we are able to detect. It is neces- 
sary to change the detection level from 5a to 3cr in order to 
observe spurious sources and make the following analysis. 
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Figure 3. Values of Iq (flux at 44 GHz) recovered by means of the MMF (asterisks) and the MF (circles). The line indicates the ideal 
recovering. The circles corresponding to the MF are slightly displaced in the horizontal axis in order to distinguish the results. 



In Figure S] we observe the number of real sources that 
both methods are capable to detect, whose intrinsic fluxes 
are higher than the corresponding value in the horizontal 
axis. As we can observe at 44 GHz, MMF detects a higher 
number of real sources for fluxes below ~ 0.4 — 0.5 Jy, being 
this difference very important at lower fluxes. Therefore, we 
obtain a clear improvement using the MMF with respect to 
the traditional matched filter at low fluxes. At 100 GHz, we 
observe a similar behaviour, but in this case the differences 
between the MF and the MMF start at ~ 0.2 Jy. If we 
observe the Figure [TJ we notice that the number of sources 
detected with the MF is higher at 100 GHz than at 44 GHz 
for values of the spectral index between and 0.5. These 
values of 7, according to the model used to simulate the 
point sources in this section, are the most frequent in the 



real sources. This gives us an idea about why the detection 
level of the MF is higher at 100 GHz. 

In Figure Owe compare the reliability of both methods 
at 44 and 100 GHz. Reliability above a certain recovered 
flux is defined as r = Nd/(Nd + N s ), where Nd is the num- 
ber of real sources above that flux, and N s is the number of 
spurious sources above the same flux. At 44 GHz we reach 
a ~ 100% of reliability at fluxes of ~ 0.3 Jy. However, the 
MF at this frequency reaches this level of reliability when 
the sources have fluxes of ~ 0.9— 1.2 Jy. At 100 GHz we ob- 
tain better levels of reliability. For example, with the MMF 
we have at 0.1 Jy more than 95% of reliability, and the MF 
reaches these values for fluxes of ~ 0.3 Jy. According to 
the expression of the reliability, this number gives us the 
percentage of real sources over the total number of sources 
detected after filtering. Therefore, we can say that the MMF 
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Figure 4. Number of real sources recovered by the MMF (solid line) and the MF (dashed line) at 44 GHz (left panel) and 100 GHz 
(right panel) whose intrinsic fluxes are higher than the corresponding value in the x axis. 



is more reliable than the MF, specially at lower fluxes. Ad- 
ditionally, we can establish the flux for which we have the 
5% of spurious sources. Making easy calculations, we finally 
obtain that the fluxes for which we have this percentage of 
spurious detections with the MF at 44 and 100 GHz are 
~ 0.5 — 0.6 Jy and 0.25 Jy respectively. If we compare these 
values with the MMF, we obtain that the fluxes are 0.15 
Jy and < 0.1 Jy. With these numbers one can see that the 
percentage of spurious detections of the MMF is much lower 
than the percentage of the MF. 

Finally, we make an additional plot where we represent, 
for both frequencies, the number of real sources detected vs. 
the number of the spurious sources (Figure [6]). In this way, 
what we represent is the number of sources that a method 
detects given a number of spurious sources. If we compare 
both plots, we can see that the curve of the MMF is always 
above the MF. It means that, when we have a fixed number 
of spurious detections, we detect more real sources with the 
MMF. 

We have to point out that the plots that we have intro- 
duced here are not directly comparable to Figure [l] As we 
have seen in this section, there are three basic and important 
differences: 

• A different way to simulate the point sources. 

• A different level of detection (in this case, a 3a level). 

• In Figure [l] we represent the number of sources with 
the corresponding flux in the horizontal axis. In the plots 
of this section, what we represent is the number of sources 
whose fluxes are higher than the corresponding value in the 
horizontal axis. 

We can say that this new multifrequency method is better in 
the sense that the number of detected sources is higher below 
~ 0.4-0.5 Jy and ~ 0.2 at 44 and 100 GHz respectively. And 
the reliability of the MMF is higher for fluxes below ~ 0.9 
Jy and ~ 0.25 Jy at 44 GHz and 100 GHz, respectively. 



5 CONCLUSIONS 

The detection of extragalactic point sources in CMB maps 
is a challenge. One has to remove them to do a proper study 



of the cosmic radiation. In addition, it is of great interest 
to study their properties, spatial and spectral distributions, 
etc. For this reason, we need suitable tools to detect and ex- 
tract these sources. There are many filtering techniques that 
have been used in this context. In this work, we have used 
the matched filter, one of the most studied techniques, and 
we have compared it with a new multifrequency one based 
on the matched multifilter (MMF). The great difference is 
that the latter takes into account information from all the 
channels of the same sky region in a simultaneous way. In 
particular, we show an example for N — 2. 

The different tests that we have used have shown an im- 
provement in the results obtained by the MMF with respect 
to the traditional matched filter. The number of detections is 
always higher when the MMF is used. In Figure[l]we see that 
we have a high number of detections with the MMF, even 
for small values of Io- It should be studied in more detail, 
but it is easy to see that one could detect and characterize 
point sources with low fluxes for 7 < 0. For this reason, this 
tool is a powerful technique to detect faint sources in CMB 
maps. 

Another important aspect is to give a good estimation of 
the quantities that we have chosen to determine the sources, 
basically the spectral index and the flux at the reference 
frequency. In both cases, we can see that the MMF improves 
the results obtained with the matched filter: the values are 
close to the input values with smaller error bars (with one 
exception, the determination of the spectral index for Io < 
0.5 Jy at positive values of the input 7). This is a significant 
fact in order to be able to detect and study properly these 
kind of sources. 

Additionally, we have made a set of more realistic sim- 
ulations in order to study and compare both filters in the 
sense of the spurious sources. We have also changed the 
threshold detection from 5a to 3a to find more spurious 
sources and make a more complete statistical analysis. First 
of all, we compare the number of real detections that we 
obtain with both techniques at 44 and 100 GHz. Compar- 
ing the plots of the Figure |4l we appreciate that, at lower 
fluxes, we detect more real sources with the matched multi- 
filter than with matched filter. This aspect is more notorious 
at 44 GHz. 
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Figure 5. Reliability versus recovered flux for the MMF (solid line) and the MF (dashed line) at 44 GHz (left panel) and 100 GHz (right 
panel). 




number of spurious defections 



number of spurious detections 



Figure 6. Number of real sources recovered by the MMF (solid line) and the MF (dashed line) at 44 GHz (left panel) and 100 GHz 
(right panel) vs. the number of spurious sources. 



One can also study the reliability of both methods. One 
can obtain a high number of real detections, and simulta- 
neously find a large number of spurious sources. Precisely, 
the reliability is a quantity that gives the number of real 
detections over the total number of sources detected. Com- 
paring the plots of the Figure 02 one can observe that the 
reliability of the matched multifilter is much higher than the 
reliability of the matched filter for low fluxes. This difference 
is particularly important at 44 GHz, where the matched fil- 
ter obtains similar values to the reliability of the matched 
multifilter only for fluxes close to 1 Jy. At 100 GHz, the 
matched filter reaches the reliability of the MMF at 3 Jy. 

The last aspect that we use to compare both methods is 
to look at the number of real sources that we have for a fixed 
number of spurious detections. The most efficient method is 
the one that has higher number of real detections for the 
same value of spurious detections. If we see the Figure [6l 
the best method is the MMF because its curves are always 
above the MF. This means that, if we take a number of 
spurious sources, the MMF recovers a larger number of real 
objects. 

Finally, we have commented at the beginning of the 



subsection 12.31 the possibility of devising a filtering method 
(the MTXF) that does not use the frequency dependence of 
the sources altogether, totally independent of the frequency 
behaviour of the sources (flat, steep or inverted). This fact 
is significant in the sense that this filtering method is a ro- 
bust technique for changes of /„. By contrast, it is necessary 
to impose the co ndition of orthonormali satio n of the matrix 
of the filters (see lHerranz fc Sand l|2008l ) and lHerranz et all 
(2009) for more details). This condition minimizes the power 
of the method. Meanwhile, the MMF is more optimal in the 
sense of the SNR (see section l2~3]l . but more complicated be- 
cause we have to maximize another set of parameters (fv). 
As one can see, the MTXF and the MMF are complemen- 
tary. 
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